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Abstract 
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of the components of the system, assuming the underlying distribution of the life¬ 
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1 Introduction 


Suppose that a series system with k repairable and identical components is under oper¬ 
ation. We assume that these k components work independently and that the lifetimes 
of all components are identically distributed. The lifetime of this system is equal to the 
minimum lifetime of its components. We also assume that the system is repaired using 
the minimal repair strategy. In a minimal repair strategy, after each failure, the corrupt 
component is immediately repaired and restored to its original condition, thus putting 
back the system into operation. Therefore, the state of the system after a repair is the 
same as it was immediately before corresponding failure. We shall assume that the time 
needed for repair and replacement is negligible. The minimal repair times possess the 
same joint distribution as upper record values from the distribution of the lifetime of the 
system, that is the distribution of record values from the distribution of the minimum of 
a sample of size k (see Barlow and Hunter, 1960). It has been verified that the sequence 
of record values, from the distribution of minimum in a sample of size k, and the sequence 
of ^-record values, from the parent distribution, are identically distributed (see Arnold et. 
ah, 1998, p. 43). 

The results of this paper focus on predicting the minimal repair times of a series system 
based on an available hybrid censored sample of its components’ lifetimes. Consider a 
sample of n units placed on a life-test at time 0. In Type-I censoring scheme, the test is 
terminated at a pre-fixed time T, while in Type-II censoring scheme, it is terminated as 
soon as a pre-determined number, r (r < n), of units fail. Under Type-I censoring scheme, 
the duration of the life-test is guaranteed, while the number of failures is random, which 
might result in a low efficiency, when the number of failures is small. In Type-II censoring 
scheme, the level of efficiency is guaranteed, as the number of failures, r, is pre-fixed, 
while the exact time of the r th failure is random, thus the duration of the experiment 
may end up being too long. The mixture of the Type-I and Type-II censoring schemes, 
called hybrid censoring scheme, proposed first by Epstein (1954), reduces the mentioned 
disadvantages. Under a Type-I hybrid censoring scheme, the experiment is terminated as 
soon as either r units fail or the time T is reached. 

Hybrid censoring has received a considerable attention in the context of reliability and 
life-testing. Many authors, including Draper and Guttman (1987), Kundu and Gupta 
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(1988), Ebrahimi (1992), Childs et al. (2003), Kundu (2007) and Kundu and Baner- 
jee (2008), have developed statistical inference based on hybrid censored sample. For a 
comprehensive review of hybrid censoring, see Balakrishnan and Kundu (2012). 

The real data used in this paper includes the number of revolutions to failure of ball 
bearings under a life test, accelerated by hybrid censoring. The ball bearings are identical, 
thus the components’ lifetimes follow the same distribution. The test is performed before 
placement of k identical ball bearings in a machine. The machine, made up of the k such 
identical components fails as soon as the first ball bearing fails, that is that the machine 
is a series system of k identical components. In the case of the failure, the physical and 
statistical (black box) minimal repair of the system are equivalent and are performed by 
minimal repair of the failed component. Our aim here is to predict the minimal repair 
times of this machine, using the information achieved from the available censored sample. 

In this paper, we assume that the underlying lifetime distribution is the two parameter 
Rayleigh distribution, with cumulative distribution function (cdf), 



( 1 ) 


where /i6l and a > 0. When /i = 0, the distribution 0 is called the scaled Rayleigh 
distribution. The corresponding probability density function (pdf) of ([!]) is 



( 2 ) 


The Rayleigh distribution is widely applied in several areas of statistics, partly be¬ 
cause of its linear and increasing failure rate, which makes it an appropriate distribution 
for modeling the lifetime distribution of components, which age rapidly with time. Several 
types of electro-vacuum devices have this feature (Polovko, 1968). The Rayleigh distribu¬ 
tion was originally introduced by Lord Rayleigh (1880, 1919) in connection with a problem 
in the field of acoustics. Wide applications of the Rayleigh distribution in lifetime analysis 
is mentioned by many authors, including Polovko (1968), Johnson et al. (1994), Dyer 
and Whisenand (1973), Gross and Clark (1976), Balakrishnan (1989) and Lawless (2003). 
The Rayleigh distribution relates to a number of well-known life distributions such as 
generalized extreme value, Weibull and Chi-square distributions (see Dey and Dey, 2014). 
There are many papers dealing with estimation and prediction under Rayleigh distribu¬ 
tion, including Howlader (1985), Howlader and Hossain (1995), Fernandez (2000), Raqab 
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and Madi (2002), Ali Mousa and Al-Sagheer (2005), Dey and Das (2007), Soliman and Al- 
Aboud (2008), Khan et al. (2010), Dey and Dey (2012) and Abdel-Hamid and Al-Hussaini 
(2014). 

Considerable work has been done on prediction of future records and order statistics, 
using parametric and nonparametric inferential methods. The following papers consider 
the Rayleigh distribution as the underlying distribution and develop predictive inferential 
methods for records and order statistics. Howlader (1985) obtained the highest posterior 
density (HPD) prediction intervals for future order statistics from an independent sample, 
based on an observed sample of order statistics. Fernandez (2000) considered the prob¬ 
lem of Bayesian prediction of a future observation based on an observed Type-II doubly 
censored sample. Raqab and Madi (2002) developed Bayesian prediction of the total time 
on a test based on doubly censored sample. Ah Mousa and Al-Sagheer (2005) considered 
Bayesian prediction of a progressive censored sample on the basis of an observed progres¬ 
sively Type-II censored sample. Soliman and Al-Aboud (2008) derived Bayesian prediction 
intervals for a future record value based on an observed sample of record values. Recently, 
Khan et al. (2010) have develpoed Bayesian inference about a future order statistic on 
the basis of a doubly censored sample. 

The observed sample and the predicted future observation might be either of the same 
type or of different types. Recently, Ahmadi and MirMostafaee (2009), Ahmadi and 
Balakrishnan (2010), Ahmadi et al. (2010) and MirMostafaee and Ahmadi (2011), have 
considered the problem of predicting future records from a sequence of observations on 
the basis of order statistics observed from another independent sample and vice versa. 
According to our best knowledge, there is a few number of similar works in the context of 
prediction of record values based on an available censored sample. These works use simple 
Type-II censored sample as the available data. The hybrid censoring is a mixture of the 
Type-I and Type-II censoring schemes. Thus the termination time of the experiment in 
hybrid censoring would decay stochastically with respect to Type-I and Type-II censoring 
schemes. 

In this paper, we obtain several Bayesian point predictors as well as Bayesian prediction 
intervals for a future repair time of a fc-component series system, applying a minimal repair 
strategy, on the basis of observed hybrid censored sample of its components lifetime, when 
the underlying distribution is assumed to be scaled or two parameter Rayleigh with cdf 
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©> and /i = 0 or p £ R, respectively. 

The rest of this paper is organized as follows. Section 2 is devoted to the description 
of the model. In Section 3, we develop the main results. Sections 4 contains a real data 
example, which illustrate the results. A simulation study is presented in Section 5 for 
comparison of the performance of the proposed predictors. 

2 The model 

Suppose that n identical components are under a life test. Let X±,..., X n denote the 
lifetimes of the experimental units and X\ n < ... < X n:n stand for the corresponding 
order statistics. Furthermore, suppose that the experiment is terminated according to a 
hybrid censoring strategy. Under a hybrid censoring scheme, the experiment is terminated 
at time Tq = minjA",.^, T}, r < n, where r and T are pre-determined values. The 
observed hybrid censored sample is therefore X = (X \. n ..., X c [, n ), where d = max{s : 
s < r, X s:n < T}. To simplify the notation, we henceforth denote the hybrid censored 
sample by (Xi,X 2 , ■■■ , X d ). 

When Xi,..., X n are independent and identically distributed (iid) with the common 
absolutely continuous cdf F and corresponding pdf /, the joint pdf of the hybrid censored 
sample is 

d 

fx u ...,x d {xi,...,x d ) = c]l f[xi) [1 - F(To)} n ~ d , (3) 

1=1 

where C is the normalizing constant. 

Suppose that the test is performed before placement of k identical components in a 
machine. Assume further that the machine, made up of the k identical components, fails 
as soon as the first component fails, that is the machine is a series system of k identical 
components. 

In the sequel, we develop Bayesian prediction of the repair times of this machine, 
repaired using a minimal repairing strategy, based on the available hybrid censored sample 
(X\. X' 2 , ■ ■ ■ , Xd). The future m th repair time of a series system, with k independent and 
identical components, the lifetimes of which have the same distribution as X\ , is denoted 

by Lm(fc) • 
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Let the sequence {Yi}X =1 be independent and identically distributed random vari¬ 
ables, independent of the sample of component lifetimes, X\.... ,X n , each of which hav¬ 
ing the same distribution as X\. Also, let Yj :n stand for the j th order statistic among 
Yi, Y n . The repair time, U m a.), is distributed as (T(m, k) — k + l) th order statistic 
among Y\, , Yr(m,k)i where T(l, k) = k and 

T(m, k) = min{j : j > T(m - 1, k), Yj > Y T(m _ i, k )-k+i:T(m-i,k)}, 

for m > 2. 

Indeed, Uu k ), U 2 (k), ■ • • is the sequence of current largest Ys yet seen (see Arnold 
et ah, 1998), that is U m n^ is identically distributed to the m th L-record value. 

The sequence {U m ^} m >i from a cdf F is identical in distribution to a record sequence 
{t/ m (i)}m>i from the cdf of the minimum in a sample of size k, F\ :k = 1 — (1 — F) k . 
Consequently, the pdf of U m ( k ) is given by (see Arnold et al., 1998) 

A'„ ( «(x;«) = log(l - F(u:e)) k r~\l - F(«; «))*-■/(«; 9). (4) 

For a comprehensive treatment and for references to the extensive literature on the 
topic of /c-record statistics, one may refer to the books of Arnold et al. (1998) and Nevzorov 
(2001). 


3 Main results 


From ([3]), the likelihood function of 9 = (//, a), under the Rayleigh distribution in ([!]), is 
given by 


L(8 |x) oc a 


—d 




_i =1 


exp 


-1 

2 a 


^2(xi - fj.) 2 + (ra - d)(T 0 - nf 
1 J I 

fi < x 1 , a > 0. (5) 

Also, from Q, the pdf of U m ^, under the Rayleigh distribution in 0 , is 


fu m(k) (u\9) = 


k m \u — n) 


2m— 1 


exp 


k(u — nY 


u > fl. 


( 6 ) 


1 ) ~~‘ jr ( 20- 

In the sequel, we provide the predictive inference, under both scaled and two parameter 
Rayleigh distributions. 
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3.1 Prediction of a repair time under the scaled Rayleigh distribution 

For the scaled Rayleigh distribution (the case p = 0), we consider the non-informative 
prior 


7Ti(cr) oc —, a > 0. 

<7 

From ([5]) and Q, the posterior density of a is given by 

-<5(x) 


^ [<Kx)] d 

7Ti (a x = ——— -z — r— exp 
iV 1 ' T(d)2 d a d + 1 


2(7 


, a > 0 , 


(7) 

( 8 ) 


where 


d 

<H X ) = + ( n ~ d)Tl (9) 

i =1 

which is the pdf of an inverted gamma distribution with parameters d and 
From ([ 8 ]) and ([ 6 ]), the predictive density of U m ^ k y given X = x, is 

/♦OO 

fu m 0 fc) H x ) = / /^ w H CT KH x ) dcJ 

d 0 

k rn u 2rn - l [5{y.)) d | ku 2 + S(x) | 

r(rn)r(d)c7 m + d + 1 2 d + m - 1 eXP l 2a J ° 



Pk(u; x) d (l - p k (u; x)) m , u > 0, 


B(d, m)u 

where B(-, •) is the complete beta function, 

1 ^ < 5 ( x ) 

and d(-) is given in ©• 

From (10), it follows that Pk{U m t k y x)|x ~ Beta(d,m), hence 

Pr{U m ( k) >z |x) = Pr(p k (U m ( k y,x) < p k (z-,x) |x) 

= I(d, m,p k (z; x)), z > 0, 

where /( 71 , 72 , x) is the incomplete beta function, 

1 


7(71,72,®) = 


>-/* 


^( 71 , 72 ) 




( 10 ) 


( 11 ) 


( 12 ) 
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Remark 1 The prior used in ([7]) is a special case of the conjugate inverted gamma prior 
with pdf 

^>0, (13) 

where a and b are positive hyper-parameters. The prior 7 r(cr) tends to 717 (<r) as (a, b) —> 
(0,0). Since then the variance of o tends to infinity, we shall call 717 (cr) the noninformative 
conjugate prior for a. 


7r <7 


oc exp 


3.1.1 Interval prediction of a repair time 


A two-sided equi-tailed 100(1 — a)% Bayesian prediction interval (PI) for U m r k \, m > 1, 
is obtained from ( 12 ), as the interval (L(x), C/(x)), for which 

(X (X 

p r( U m(k) > T(x)|x) = 1 - - and Pr(U m ( k) > U(x) |x) = -. 

So, the interval (L(x),f7(x)) is 


[1 - /3f (d, m)]d(x) / [1 - /3i_a (d, m)]5(x) 


kfiz (d, m) 


A:/3i_|(d, m) 


(14) 


where /3 7 (ni, 712 ) denotes the upper 7 th quantile of the beta distribution with parameters 
n 1 and n 2 , i.e., P(T > / 3 7 (ni,n 2 )) = 7 , with T ~ Beta(n\,n 2 )- 

The highest posterior density prediction interval (HPD PI) is an interval, the posterior 
pdf for every point inside which is greater than that for every point outside of which. 
A HPD PI includes the more probable values of the parameter and excludes the less 
probable ones. Since the posterior pdf f? T (wlx) is unimodal and »fc(u;x) is decreasing 

^m(k) x v 7 

in u. the HPD PI, (uq, W 2 ), for U m ^.\ given X = x, with coverage probability 1 — ct, is the 
simultaneous solution of 


1 

B(d, m) 


rPk(w i;x) 
'Pk(w 2;x) 


t d -\l-t) m ~ 1 dt = 1-a, 


(15) 


and 


W\ 

W2 


2m— 1 


/(5(x) + kw\ \ d+m 
\5(x.)+kwlJ 


( 16 ) 









Remark 2 For m = 1, we have U j= Yi : fc, that is the lifetime of the series system, and 
the prediction interval (14) is simplified to 



Also, the HPD PI, (ici, W 2 ), for given X = x, with coverage probability 1 — a, is the 
simultaneous solution of 

\Pk(w i; x)] d - \pk{w 2 ; x)] J = 1 - a 

and 

^i[Pfc(^i;x )] d+1 = w 2 \p k (w 2 -,x)] d+1 . 


3.1.2 Point prediction of a repair time 


Using (10) and under the squared error loss (SEL) function, the Bayes point predictor of 


Um(k) |x) 

roo 2 


f 


B(d, m 


-Pk(w, x ) d ( 1 - Pk(u ; x)) m du 


5(x) 

k z d ~l(l- z) m ~hz 


J o B(d,m 
B(d-l,m+±) J5(X) 

B(d,m) V k ' ( ’ 

Similarly, it can be verified that, under the absolute error loss (AEL) function and zero-one 
loss function, the Bayes point predictors of U m ^) are 


U* m{k) = Med(t/ m(fc) |x) = 


(1 — Med(d, m))[<5(X)] 


fcMed(d, m) 


and 


U m(k) = Mod(V : „ w |x) = '/ (2 ” (M 1 f 1 ( ) X>l 


respectively, where Med(d, m) denotes the median of Beta distribution with parameters d 
and m. 
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3.2 Results for the two parameter Rayleigh distribution 

To facilitate the Bayesian approach under the two parameter Rayleigh distribution, we 
assume independent prior distributions for the model parameters, that is 


^ r(/i,cr) oc 7ri(cr)7r 2 (», 


(19) 


where tt\ (cr) is the non-informative conjugate prior in ([7]) and 7r 2 (/r) is a normal density 
with mean £ and variance l/2r, that is 


7t 2 (/x) oc exp {— t(h — £) 2 } , |i£l, £ 6 M, r > 0. 


Therefore, the joint prior density is 

1 
cr 


vr 


(H, cr) oc — exp {—r(/x — £) 2 ) , cr > 0, /jGE, (gR, r > 0, 
a y ’ 

and the posterior density function of g and cr, given X = x, is obtained as 
Ai(x)exp{—T(/r-£) 2 } r d 


7T*(^,cr|x) = 


T(d)2 d a d+1 


Y[(xi - n) 


_i=l 


ex P { ^[<5*Ol x )] 


cr > 0, fi < X i, 


where 


<**(mI x ) = ~ df + {n- d)(T 0 - gf 


2—1 


and 


^i( x ) 1 = 


f Xl exp{-r(/r-0 2 }nti(^-^) 


dg. 


[<5*(//|x)] d 

The predictive posterior density of U m ^, given x, is obtained as follows 

f*min(2i,ici) roc 


l U 


m(k) 


/ min^tx,£1; roo 

/ cr)7T*(At, cr| x )dcrd/x 

-00 J 0 


2 Ai (x) /• min («^i) g (t,u,m-x,r]) 


u — t 


df, u E 


B(d, m) 7„ oc 
where r] = (r, £, k, d) and 

k m exp {—r(f - £) 2 } (u ~ t) 2m Y\i=i( x i ~ t) 


g(t,u,m-,x,rj) = 


[k(u — t) 2 + 5*(t |x)] d+m 


( 20 ) 


( 21 ) 


( 22 ) 


(23) 


(24) 


(25) 


(26) 
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The predictive posterior survival function of U m n^ is, for z > x i, 


/ , X\ /* oo /*oo 

/ / /t^ooM/^^O^dudad/i 

-oo J 0 J z 
m— 1 

= Ar( 

j=o 


, (d + j) f xi , 

L s{t ’ z ' r '*’ ri)dt ' 


( 27 ) 


where g(t,z,j\-x.,rj) is given in (26); and for z < x i 


S u mW ( z W = ^ x ) 




+ g(t, t, 0 ; x, 77 ) df 

J z 

wherein 0 ° is defined to be 1 and 5 is given in (26). 


(28) 


3.2.1 Interval prediction of a repair time 

The equi-tailed 100(1 — a)% Bayesian prediction interval for U m ^\ can be obtained nu¬ 


merically, using (27) and (28). 

Since hy (u|x) is unimodal, then the HPD PI, ( 101 , 102 ), for U m ^, with coverage 
probability 1 — a, satisfies 


rw 2 


f W 1 


hu , s (u\x.)du = 1 — a 

^rn{k) x ' 


and 


h U mW ( w i\*) = h U mW Mx), 


where hr r (ulx) is given in (25). 

^rn(k) v 7 - 


3.2.2 Point prediction of a repair time 

Under the SEL function, the point predictor of U m ^) is given by 

Um.(k) E(Um(k) 1^0 

T(m + l/2)T(d — 1/2) f Xl 


= Ai(x) 
rx 1 


T(m)T(d)Vk 

/ Xi 

t 3 (t,t, 0 ;x,T 7 ) df 

-00 


[ V^*(i\x)g(t, t, 0 ; X, 77 ) df 

J — OO 


( 29 ) 


where g is given in (26). 
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For the AEL function, U* m ^ can be obtained as 


U* 


= H*u 

rr 


™( k ) — ^m(fc) l 2 

Also, the point predictor U m ( k )i under the zero-one loss function is the unique mode of 


the pdf hrr (re lx) in (25). 

u m(k) x 7 - 


3.3 Model checking 


The assumption of independence of 7 Ti(<t) and vr 2 (^) in (19) may affect on the perfor¬ 


mance of the predictors. Furthermore, as r tends to zero, the prior (m) tends to a 
noninformative prior, while the mean square error of prediction, 

^(^ro(fc) ^rn(k )) i 

increases. Hence, choosing a suitable value for r is an important issue. It should be kept 


in mind that U ., 


m(k) 


minimizes 


E[(g(x) - G m(fc) ) 5 


x 


over all functions g, while E (U m ^) ~ ^m(fc)) 2 can be greater or less than — 

U m (k )) 2 = 2Var([/ m (fc)), in which U rn ^ and U' m t k ) are iid random variables from the 
distribution Q. Hence,for the prior distribution in (19) to be sufficiently low informative, 
and for the mean square error of prediction not to be very large, one may choose r such 
that ^ 

~ U m ( k )) 2 

takes its largest value less than or equal to 1. 


To check the suitableness of the prior distribution in (19) and to choose a suitable 


value for r, we perform a simulation study as follows. This study is based on the general 
method for model checking, described by Gelman et al. (2004). 

The algorithm for model checking is as follows: 

Algorithm: 

(i) Generate aq, ■ ■ ■ ,x n independently from Rayleigh distribution and extract the hybrid 
censored sample (aq, • • • , Xd) from aq, • ■ ■ ,x n . 
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(ii) Generate u m ^) and u', k , independently from 


(iii) Predict u* m ^ using p9[ ). 

(iv) Replicate (i)-(iii), N = 10000 times, independently, to obtain samples u^} k y 


and u'^ k y ..., u'^ k j as well as predicted values u*' m ^ k y ■ ■ ■, ur m ( k y 


1 .( 1 ) 




IN) 

^m^k) 


(v) Compute 


and 


N 


SS ' = N^ 




1.(0 


^ m(k) ^ m(k)) 


i= 1 


N 


SS 2 = ^ - w 


*m(k) 


'(*) l 2 

m(k)' ’ 


i =1 


(vi) Set £ = 0 and choose l* from the set {— 2 , — 1 ,..., 2 }, such that, for r* = 0.5 x 10 
D\ = SS 1 /SS 2 takes its largest value less than or equal to 1 . 


(vii) For £ = 0 and r = r*, compute 


£>2 = 


at 


'm(fc) 




m(k) 


-M) 


and 


iESAMd2i». 


H 3 = 


i(fc) 


l m(k) 


lyW (».*(*) \2 _ / J_ 

AT l at 2Gj=i 


rt' 


m(k) 


2 ‘ 


Table [l] presents the values of r*, Hi, H 2 and H 3 , for m = 3, k = 2, yu = 0, a = 1, 
£ = 0 and r = r*. Also, empirical cdfs of the simulated values, u^/ k y ■ ■ ■ i u m( k y and the 
predicted values, u*^ k y... ,u*^ k y are shown in Figure |TJ for m = 3, k = 2 , fi = 0 , 
a = 1, £ = 0, r = r* and different values of n and (r, T). As one can see from Table [I] 
and Figure [lj for a sufficiently low informative prior, the predicted values have smaller 
variance and mean square error of prediction than the simulated values. The ignorable 


bias of the predicted values decreases as n and/or r get large. Hence, the prior (19) results 
in efficient predictors. 
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Figure 1: Empirical cdfs of the simulated values, ■ ■ ■ ■> u m(k) ( s °dd line), and the 

predicted values, u*^ k y ■ ■ ■, u *^ k ) (dashed line), for m = 3, k = 2, ii = 0, a = 1, £ = 0, 
t = t* and different values of n and (r, T). ^ 































Table 1: The values of r*, D\, D -2 and D 3 , for m = 3, k = 2, /i = 0, a = 1, £ = 0 and 


T = T*. 


n 

(r,T) 

T* 

Di 

D-2 

d 3 

10 

(5,2) 

0.5 

0.9426 

0.9330 

1.2006 


(5,2.5) 

0.5 

0.9417 

0.9348 

1.2513 


(8,2) 

0.005 

0.8605 

0.9441 

1.5197 


(8,2.5) 

0.005 

0.7295 

0.9508 

2.2776 

20 

(10,2) 

0.05 

0.7160 

0.9577 

2.5605 


(10,2.5) 

0.005 

0.7143 

0.9561 

2.4767 


(16,2) 

0.005 

0.6188 

0.9750 

4.3192 


(16,2.5) 

0.005 

0.6125 

0.9806 

5.1078 

30 

(15,2) 

0.05 

0.6302 

0.9740 

4.2603 


(15,2.5) 

0.005 

0.6262 

0.9758 

4.1887 


(24,2) 

0.005 

0.5595 

0.9852 

6.9719 


(24,2.5) 

0.005 

0.5714 

0.9858 

7.8098 


4 Illustrative example: groove ball bearings data set 

In this section, we illustrate the proposed procedures in the previous section, using the 
following real data example. 

Example 4.1 Consider the endurance test of deep groove ball bearings, discussed by 
Leiblein and Zelen (1956), which includes the number of revolutions (in hundreds of mil¬ 
lions) to failure, for each of n = 23 ball bearings, as follows: 

0.1788 0.2892 0.3300 0.4152 0.4212 0.4560 0.4880 0.5184 0.5196 0.5412 
0.5556 0.6780 0.6864 0.6864 0.6888 0.8412 0.9312 0.9864 1.0512 1.0584 
1.2792 1.2804 1.7340. 

The adequacy of the fitness of the two parameter Rayleigh distribution with fi = 0.1788 
and <7 = 0.2149 to the data is tested using the Kolmogorov-Smirnov (K-S) test. The 
value of K-S test statistic is obtained as D = 0.1982 with a corresponding j>value=0.3269. 
Hence, the two parameter Rayleigh distribution fits the data quite well. 

Consider the following two sampling schemes: 
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• Scheme 1: r = 20 and T = 1.25. 

We have, d = r = 20, To = min{x r:n ,T} = 1.0584, and the hybrid censored sample 

is 

0.1788 0.2892 0.3300 0.4152 0.4212 0.4560 0.4880 0.5184 0.5196 0.5412 

0.5556 0.6780 0.6864 0.6864 0.6888 0.8412 0.9312 0.9864 1.0512 1.0584. 

• Scheme 2: r = 20 and T = 1. 

We have, in this case, d = 18, To = 1 and the hybrid censored sample is 

0.1788 0.2892 0.3300 0.4152 0.4212 0.4560 0.4880 0.5184 0.5196 

0.5412 0.5556 0.6780 0.6864 0.6864 0.6888 0.8412 0.9312 0.9864. 

To obtain stable Bayesian predictors, we need to choose suitable values for the hyper¬ 
parameters £ and r, which ensures that the prior distribution is sufficiently low-informative. 
This is performed through a sensitivity analysis in the next subsection. 

4.1 Sensitivity analysis 

In order to check the effect of hyper-parameters £ and r on Bayesian predictors, we consider 
a pilot run on the prediction procedure, for various values of £ and r. As r —> 0, the 
variance of the normal prior for the parameter /i tends to infinity, which is the non- 
informative case. As the variance of the normal prior increases, the effect of the hyper¬ 
parameter £ on the predictors decreases. Therefore, many authors perform the sensitivity 
analysis on the hyper-parameter r, assuming £ to be fixed as £ = 0. So, we take £ = 0 and 
focus on a sequence of values of r, tending to 0. Figure [ 2 ] shows the plot of U m ^) versus 
values of l, for r = 0.5 X 10 _i , and different values of m and k , when £ = 0. Figure [ 2 ] 
shows that for r > 0.5, the values of U m ^ tend to be very close to each other. Therefore, 
r = 0.5 and £ = 0 are suitable hyper-parameters for obtaining stable Bayesian predictors. 

4.2 Bayesian prediction 

Here, for the sake of comparison of the performance of the predictors, we take k = 1,2,3 
and m = 1,..., 4. 
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Table 2: 95% Bayesian prediction intervals and point predictors for U m i k \, based on the 


hybrid censored sample in Example 4.1 


Scheme 1 



m 

Equi-tailed PI [width] 

HPD PI [width] 

^ra(fc) 

U* 

m(K) 



i 

(0.1575,1.5235) [1.3660] 

(0.1004,1.4219) [1.3215] 

0.7111 

0.6663 

0.5692 


2 

(0.4070,1.8904) [1.4834] 

(0.3468,1.7975) [1.4507] 

1.0402 

1.0027 

0.9302 

k = 1 

3 

(0.6175,2.1717) [1.5542] 

(0.5567,2.0811) [1.5244] 

1.2870 

1.2498 

1.1796 


4 

(0.7965,2.4113) [1.6148] 

(0.7337,2.3184) [1.5847] 

1.4928 

1.4542 

1.3825 


1 

(0.1105,1.0894) [0.9789] 

(0.0736,1.0329) [0.9593] 

0.5183 

0.4895 

0.4260 


2 

(0.2966,1.3448) [1.0482] 

(0.2580,1.2885) [1.0305] 

0.7510 

0.7268 

0.6798 

k = 2 

3 

(0.4495,1.5408) [1.0913] 

(0.4099,1.4839) [1.0740] 

0.9255 

0.9014 

0.8556 


4 

(0.5790,1.7078) [1.1288] 

(0.5377,1.6485) [1.1108] 

1.0710 

1.0459 

0.9989 


1 

(0.0851,0.8988) [0.8137] 

(0.0584,0.8616) [0.8032] 

0.4329 

0.4114 

0.3632 


2 

(0.2445,1.1051) [0.8606] 

(0.2160,1.0656) [0.8496] 

0.6229 

0.6048 

0.5693 

k = 3 

3 

(0.3724,1.2633) [0.8909] 

(0.3424,1.2219) [0.8795] 

0.7654 

0.7472 

0.7124 


4 

(0.4801,1.3982) [0.9181] 

(0.4485,1.3542) [0.9057] 

0.8842 

0.8650 

0.8291 

Scheme 2 


m 

Equi-tailed PI [width] 

HPD PI [width] 

^ra(fc) 

77* 

m(k) 



i 

(0.1487,1.5871) [1.4384] 

(0.0880,1.4792) [1.3912] 

0.7291 

0.6811 

0.5785 


2 

(0.4101,1.9789) [1.5688] 

(0.3444,1.8754) [1.5310] 

1.0737 

1.0320 

0.9526 

k = 1 

3 

(0.6290,2.2807) [1.6517] 

(0.5612,2.1767) [1.6155] 

1.3322 

1.2898 

1.2112 


4 

(0.8143,2.5383) [1.7240] 

(0.7435,2.4299) [1.6864] 

1.5476 

1.5031 

1.4215 


1 

(0.0980,1.1290) [1.0310] 

(0.0592,1.0700) [1.0108] 

0.5272 

0.4967 

0.4302 


2 

(0.2947,1.4011) [1.1064] 

(0.2527,1.3390) [1.0863] 

0.7709 

0.7442 

0.6931 

k = 2 

3 

(0.4543,1.6108) [1.1565] 

(0.4102,1.5461) [1.1359] 

0.9536 

0.9263 

0.8754 


4 

(0.5888,1.7901) [1.2013] 

(0.5422,1.7212) [1.1790] 

1.1060 

1.0771 

1.0238 


1 

(0.0706,0.9282) [0.8576] 

(0.0427,0.8900) [0.8473] 

0.4377 

0.4153 

0.3653 


2 

(0.2398,1.1474) [0.9076] 

(0.2090,1.1045) [0.8955] 

0.6367 

0.6169 

0.5786 

k = 3 

3 

(0.3737,1.3164) [0.9427] 

(0.3405,1.2697) [0.9292] 

0.7860 

0.7655 

0.7269 


4 

(0.4860,1.4609) [0.9749] 

(0.4504,1.4101) [0.9597] 

0.9103 

0.8884 

0.8479 
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Figure 3; Predictive density of m th repair time of a series machine with k ball bearings (in hundreds of millions 
revolutions) given the censored sample of Example |4.l| (Scheme 1) for k = 1,2,3 (rows from top to bottom) and 
m = 1, 2, 3 (columns from left to right). 
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We shall construct equi-tailed 95% Pis, as well as 95% HPD Pis, for the first four repair 
times (in hundreds of millions revolutions to failure), based on a minimal repair strategy, 
of a series system with k = 1, 2, 3 ball bearings, i.e. for U m ^, m = 1,..., 4 and k = 1,2,3. 
The intervals and the corresponding widths, as well as the point predictors U m ^, 17* ^ 
and are calculated and presented in Table [IJ for each of the two sampling schemes. 

For example, for k = 3 and m = 4, the point predictor 1 / 4 ( 3 ) = 0.9103 means that, when 
a minimal repair strategy is used, on average a series system with 3 ball bearings fails 
and needs to be repaired for the fourth time after 91.03 million revolutions. One can see 
from Table [2] that HPD Pis are more precise than their corresponding equi-tailed Pis. We 
also provide the predictive density plots of U m ^, for k = 1,2,3 and m = 1, 2,3, based on 
Scheme 1 in Figure [3} 


5 Simulation study 


In this section, we wish to compare the performance of the interval and point predictors 
through a simulation study. As the Rayleigh distribution belongs to location-scale family, 
it would be reasonable to take /r = 0 and a = 1 in the simulation. To examine the 
effect of hyper-parameters on the predictors, the informative prior is considered by setting 
£ = —1,1 and t = 0.25, 5. Moreover, we consider the low-informative case by setting £ = 0 


and t = 0.005 in (21). We take n = 20, r = 17 T = 2, k = 1,2,3 and m = 1,2,3. The 
following algorithm is used to perform the simulation: 


1. Generate aq, ■ ■ ■ ,x n independently from Rayleigh distribution and extract the hybrid 
censored sample (aq, • • • , xj) from aq, ■ ■ • ,x n . 

2 . Generate the repair time u m (*,) from ©• 

3. Obtain 95% Pis as well as the point predictors u m (fc), and based 

on the values of (aq, • • • , Xd)- 

4. Repeat Steps 1-3 N = 10000 times, to obtain and 

(0 ; ^ f 1 • • • ; ^• 
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5. Calculate the estimated risk (ER) of the point predictors as follows 


N 


ER( u m{k)) jy ^ ' [^m(fe) (*) (*)] ) 


1 


i =1 
TV 


mK { k)) = X X K*)(o - “».<*)« 


and 


where 


2—1 

AT 


ER( u m(k)) y ^ M m(k )(*)) 


i=l 


5(f) = 


1, t = 0 

0, t^O. 


6. Calculate the average width (AW) and coverage probability (CP) of the Pis. 

In order to compare the Bayesian predicts with the classical (frequentist) ones, we also 
consider a classical method, due to Wald (1942), to obtain prediction intervals. To this end, 
we have to plug in the maximum likelihood estimators (MLEs) of the parameters /r and a 
based on the hybrid censored sample in the conditional (predictive) density /t/ m(fe) |x(^| x ) 
to estimate it and then use it to obtain the PI for U m ^y Using the likelihood function 
(§, the maximum likelihood estimates /j and a are obtained by solving the likelihood 
equations 


X 

2—1 


-1 2 1 
-4- 

Xi — n a 


^2(xi - n) + (n - d)(T 0 - fi) 


_ i— 1 


= o 


and 


+ 2 = °> 


where <5*(^|x) is given in (23). 

Since the hybrid censored sample X and U m ^) are independent, we have fu, n{k) \x.( u l x ) = 
f Um(k) {u) and thus the estimated predictive density is obtained as 


fu m(k) (u) = 


k m (u — fi) 


2m—1 


T(m)a m 2 


m —1 


exp 


k(u — jiy 
2d 


u > fi. 


(30) 
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Letting, U* i(k) = k(U m(fc) - fi) 2 /a, we have 

^ u m—l e ~u/2 

f u ™(k) l x ( n l x ) — r(m) 2 m_1 ’ 

that is the pdf of chi-square distribution with 2 m degrees of freedom. Thus, a two-sided 
equi-tailed 100(1 — a)% Wald’s prediction interval for U m ^, given the hybrid censored 
sample is 



where y 2 ( 2 rn) denotes the upper 7 th quantile of the chi-square distribution with 2 m de¬ 
grees of freedom, i.e., P(T > % 2 (2 m)) = 7 , with T ~ y 2 (2m). 

The simulated AWs, as well as CPs, of the Wald’s Pis for U m ^ are given in Table 
8, for k = 1,2,3 and m = 1,2,3. Also, the simulated AWs and CPs of the equi-tailed 
Bayesian and HPD Pis for U m n^ are given in Tables 3-7, for different values of the hyper¬ 
parameters £ and r, k = 1, 2,3 and m = 1, 2, 3. As one can see from Tables 3-7, the ERs 
of the point predicts and AWs of the Pis are increasing with respect to m, when other 
parameters are kept fixed, while they are decreasing with respect to k. Also we observe 
that although HPD Pis have smaller AWs in comparison with equi-tailed Pis, their CPs 
are very close and even in most cases larger than the CPs of equi-tailed Pis, revealing the 
superiority of HPD Pis. For £ = — 1 and r = 5, the largest CPs and AWs are obtained, 
while for £ = 1 and r = 5, we obtain the smallest CPs and AWs. However, the ERs of all 
cases are almost equal. From Table 8 , we see that the Wald’s Pis have smaller AWs and 
CPs in comparison with their corresponding Bayesian Pis. 

6 Concluding Remarks 

In this paper, we have studied the marginal Bayesian prediction of repair times of a series 
system based on a minimal repair strategy, using the information contained in an observed 
hybrid censored sample of the lifetimes of the components of the system. The results can 
be extended to non-parametric prediction procedures and also the joint prediction of repair 
times. The problem of predicting repair times of other types of coherent systems, such 
as parallel systems and fc-out-of-n systems, on the basis of a censored sample of their 
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Table 3: The simulated AWs and CPs of the Bayesian Pis for U m ^, for £ = — 1 and r = 5. 




m 

Equi-tailed 

AW 

PI 

CP 

HPD 

AW 

PI 

CP 

ER(U m(k) ) 


ER(U m(k) ) 



1 

3.5841 

0.9930 

3.4763 

0.9875 

0.4502 

0.5360 

1.0000 

k = 

1 

2 

3.8944 

0.9890 

3.7825 

0.9889 

0.5946 

0.5970 

1.0000 



3 

4.1232 

0.9790 

3.9736 

0.9886 

0.7854 

0.6697 

1.0000 



1 

2.5732 

0.9858 

2.5323 

0.9795 

0.2396 

0.3894 

1.0000 

k = 

2 

2 

2.7572 

0.9920 

2.6778 

0.9876 

0.2736 

0.4129 

1.0000 



3 

2.8789 

0.9482 

2.8524 

0.9741 

0.5328 

0.5600 

1.0000 



1 

2.1586 

0.9831 

2.1301 

0.9750 

0.1800 

0.3407 

1.0000 

k = 

3 

2 

2.2668 

0.9909 

2.2410 

0.9858 

0.1807 

0.3399 

1.0000 



3 

2.3508 

0.9913 

2.3073 

0.9879 

0.2029 

0.3542 

1.0000 


Table 4: The simulated AWs and CPs of the Bayesian Pis for U m n^ 

, for £ = 

1 and r = 5 




Equi-tailed 

PI 

HPD 

PI 






m 

AW 

CP 

AW 

CP 

ER(.U m{k) ) 


ER(U m(k) ) 



1 

2.4006 

0.9036 

2.2932 

0.9108 

0.4687 

0.5423 

1.0000 

k = 

1 

2 

2.6305 

0.9204 

2.5956 

0.9278 

0.5263 

0.5771 

1.0000 



3 

2.7528 

0.9222 

2.6601 

0.9195 

0.5698 

0.6023 

1.0000 



1 

1.7006 

0.8899 

1.6687 

0.9094 

0.2369 

0.3865 

1.0000 

k = 

2 

2 

1.8560 

0.9094 

1.8047 

0.9172 

0.2686 

0.4127 

1.0000 



3 

1.9450 

0.9223 

1.9250 

0.9578 

0.2848 

0.4241 

1.0000 



1 

1.4001 

0.8804 

1.3587 

0.9010 

0.1622 

0.3195 

1.0000 

k = 

3 

2 

1.5176 

0.9084 

1.5041 

0.9218 

0.1772 

0.3344 

1.0000 



3 

1.5873 

0.9187 

1.5257 

0.9249 

0.1926 

0.3471 

1.0000 
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Table 5: The simulated AWs and CPs of the Bayesian Pis for U . m (^, for £ = — 1 and 


r = 0.25. 




m 

Equi-tailed 

AW 

PI 

CP 

HPD 

AW 

PI 

CP 

ER(U m ( fc) ) 


ER(U m(k) ) 



1 

2.8046 

0.9634 

2.7415 

0.9654 

0.4619 

0.5429 

1.0000 

k = 

1 

2 

3.0324 

0.9507 

2.9718 

0.9561 

0.5473 

0.5845 

1.0000 



3 

3.2103 

0.9491 

3.1510 

0.9593 

0.6153 

0.6157 

1.0000 



1 

1.9971 

0.9570 

1.9859 

0.9582 

0.2256 

0.3784 

1.0000 

k = 

2 

2 

2.1387 

0.9523 

2.0747 

0.9524 

0.2657 

0.4101 

1.0000 



3 

2.2313 

0.9519 

2.1790 

0.9590 

0.2935 

0.4252 

1.0000 



1 

1.6724 

0.9541 

1.6606 

0.9575 

0.1598 

0.3170 

1.0000 

k = 

3 

2 

1.7622 

0.9571 

1.7524 

0.9585 

0.1764 

0.3341 

1.0000 



3 

1.8334 

0.9489 

1.8074 

0.9541 

0.1997 

0.3532 

1.0000 


Table 6: 

The simulated AWs and CPs of the Bayesian Pis for U m ^, for £ = 1 and 

r = 0.25. 











Equi-tailed 

PI 

HPD 

PI 





m 

AW 

CP 

AW 

CP 

ER(U m(k) ) 


ER(U m(k) ) 


1 

2.7305 

0.9587 

2.5025 

0.9584 

0.4628 

0.5431 

1.0000 

k = 1 

2 

2.9582 

0.9428 

2.8333 

0.9485 

0.5474 

0.5795 

1.0000 


3 

3.1104 

0.9498 

3.0390 

0.9563 

0.6083 

0.6163 

1.0000 


1 

1.9378 

0.9467 

1.9217 

0.9516 

0.2313 

0.3827 

1.0000 

k = 2 

2 

2.0774 

0.9492 

2.0129 

0.9474 

0.2632 

0.4070 

1.0000 


3 

2.1777 

0.9465 

2.1290 

0.9525 

0.2942 

0.4286 

1.0000 


1 

1.6144 

0.9482 

1.5954 

0.9483 

0.1582 

0.3173 

1.0000 

00 

II 

2 

1.7069 

0.9466 

1.6974 

0.9527 

0.1754 

0.3332 

1.0000 


3 

1.8085 

0.9633 

1.7874 

0.9678 

0.1986 

0.3577 

1.0000 
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Table 7: The simulated AWs and CPs of the Bayesian Pis for U m ^\, for £ = 0 and 
r = 0.005._ 




m 

Equi-tailed 

AW 

PI 

CP 

HPD 

AW 

PI 

CP 

ER(U m(k) ) 


ER(U m(k) ) 



1 

2.7430 

0.9471 

2.6487 

0.9486 

0.4661 

0.5379 

1.0000 

k = 

1 

2 

2.9953 

0.9489 

2.9414 

0.9539 

0.5421 

0.5791 

0.9999 



3 

3.1791 

0.9538 

3.1284 

0.9594 

0.6086 

0.6118 

0.9997 



1 

1.9735 

0.9494 

1.9649 

0.9548 

0.2317 

0.3832 

1.0000 

k = 

2 

2 

2.1178 

0.9452 

2.0502 

0.9470 

0.2696 

0.4107 

1.0000 



3 

2.2232 

0.9513 

2.1658 

0.9559 

0.2960 

0.4293 

1.0000 



1 

1.6489 

0.9493 

1.6334 

0.9525 

0.1590 

0.3173 

1.0000 

k = 

3 

2 

1.7390 

0.9526 

1.7301 

0.9571 

0.1733 

0.3313 

1.0000 



3 

1.8103 

0.9503 

1.7897 

0.9538 

0.1974 

0.3482 

1.0000 


Table 8: The simulated AWs and CPs of the Wald’s Pis for U m ^y 




m 

AW 

CP 



1 

2.3382 

0.9059 

k = 

1 

2 

2.4778 

0.9091 



3 

2.5293 

0.8943 



1 

1.6496 

0.8964 

k = 

2 

2 

1.7597 

0.9027 



3 

1.7807 

0.8978 



1 

1.3463 

0.8801 

k = 

3 

2 

1.4294 

0.8994 



3 

1.4734 

0.9022 
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components is of interest. Work on these problems is under progress and we hope to 
report the new results in the near future. 
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